In [1]:
%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
import netCDF4 as nc
import climlab
In [2]:
# Test in a 2-layer atmosphere
col = climlab.GreyRadiationModel(num_lev=2)
print col
In [3]:
col.subprocess
Out[3]:
In [4]:
col.state
Out[4]:
In [5]:
col.Ts
Out[5]:
In [6]:
col.Ts[:] = 288.
col.Tatm[:] = np.array([275., 230.])
col.state
Out[6]:
In [7]:
LW = col.subprocess['LW']
print LW
In [8]:
LW.absorptivity
Out[8]:
In [9]:
LW.absorptivity = 0.58377
LW.absorptivity
Out[9]:
In [10]:
col.diagnostics
Out[10]:
In [11]:
col.compute_diagnostics()
col.diagnostics
Out[11]:
In [12]:
col.diagnostics['OLR']
Out[12]:
In [13]:
col.state
Out[13]:
In [14]:
col.step_forward()
In [15]:
col.state
Out[15]:
In [16]:
# integrate out to radiative equilibrium
col.integrate_years(2.)
In [17]:
col.diagnostics['ASR'] - col.diagnostics['OLR']
Out[17]:
In [18]:
# Compare these temperatures against our analytical solutions for radiative equilibrium
col.state
Out[18]:
In [19]:
ncep_url = "http://www.esrl.noaa.gov/psd/thredds/dodsC/Datasets/ncep.reanalysis.derived/"
ncep_air = nc.Dataset( ncep_url + "pressure/air.mon.1981-2010.ltm.nc" )
level = ncep_air.variables['level'][:]
lat = ncep_air.variables['lat'][:]
zstar = np.log(level/1000)
In [20]:
Tzon = np.mean(ncep_air.variables['air'][:],axis=(0,3))
Tglobal = np.average( Tzon , weights=np.cos(np.deg2rad(lat)), axis=1) + climlab.constants.tempCtoK
In [21]:
fig = plt.figure( figsize=(8,6) )
ax = fig.add_subplot(111)
ax.plot( Tglobal , level )
ax.invert_yaxis()
ax.set_xlabel('Temperature (K)', fontsize=16)
ax.set_ylabel('Pressure (hPa)', fontsize=16 )
ax.set_title('Global, annual mean sounding from NCEP Reanalysis', fontsize = 24)
ax.grid()
In [22]:
# initialize a grey radiation model with 30 levels
col = climlab.GreyRadiationModel()
print col
In [23]:
# interpolate to 30 evenly spaced pressure levels
lev = col.lev
In [24]:
Tinterp = np.flipud(np.interp(np.flipud(lev), np.flipud(level), np.flipud(Tglobal)))
Tinterp
Out[24]:
In [25]:
# Initialize model with observed temperatures
col.Ts[:] = Tglobal[0]
col.Tatm[:] = Tinterp
In [26]:
def plot_sounding(collist):
color_cycle=['r', 'g', 'b', 'y']
# col is either a column model object or a list of column model objects
if isinstance(collist, climlab.Process):
# make a list with a single item
collist = [collist]
fig = plt.figure()
ax = fig.add_subplot(111)
for i, col in enumerate(collist):
ax.plot(col.Tatm, col.lev, color=color_cycle[i])
ax.plot(col.Ts, climlab.constants.ps, 'o', markersize=12, color=color_cycle[i])
ax.invert_yaxis()
ax.set_xlabel('Temperature (K)')
ax.set_ylabel('Pressure (hPa)')
ax.grid()
return ax
In [27]:
plot_sounding(col)
Out[27]:
In [28]:
col.compute_diagnostics()
col.diagnostics['OLR']
Out[28]:
In [29]:
# Need to tune absorptivity to get OLR = 239
epsarray = np.linspace(0.01, 0.1, 100)
OLRarray = np.zeros_like(epsarray)
In [30]:
for i in range(epsarray.size):
col.subprocess['LW'].absorptivity = epsarray[i]
col.compute_diagnostics()
OLRarray[i] = col.diagnostics['OLR']
In [31]:
plt.plot(epsarray, OLRarray)
plt.grid()
In [32]:
def OLRanom(eps):
col.subprocess['LW'].absorptivity = eps
col.compute_diagnostics()
return col.diagnostics['OLR'] - 239.
In [33]:
OLRanom(0.02)
Out[33]:
In [34]:
# Use numerical root-finding to get the equilibria
from scipy.optimize import brentq
# brentq is a root-finding function
# Need to give it a function and two end-points
# It will look for a zero of the function between those end-points
eps = brentq(OLRanom, 0.01, 0.1)
print eps
In [35]:
col.subprocess['LW'].absorptivity = eps
col.subprocess['LW'].absorptivity
Out[35]:
In [36]:
col.compute_diagnostics()
col.diagnostics['OLR']
Out[36]:
In [37]:
col2 = climlab.process_like(col)
print col2
In [38]:
col2.subprocess['LW'].absorptivity *= 1.02
col2.subprocess['LW'].absorptivity
Out[38]:
In [39]:
col2.compute_diagnostics()
col2.diagnostics['OLR']
Out[39]:
In [40]:
col2.Ts - col.Ts
Out[40]:
In [41]:
col2.diagnostics['OLR'] - col.diagnostics['OLR']
Out[41]:
In [42]:
RF = -(col2.diagnostics['OLR'] - col.diagnostics['OLR'])
print 'The radiative forcing is %f W/m2.' %RF
In [43]:
re = climlab.process_like(col)
In [44]:
re.integrate_years(2.)
In [45]:
# Check for energy balance
re.diagnostics['ASR'] - re.diagnostics['OLR']
Out[45]:
In [46]:
plot_sounding([col, re])
Out[46]:
In [47]:
rce = climlab.RadiativeConvectiveModel(adj_lapse_rate=6.)
print rce
In [48]:
rce.subprocess['LW'].absorptivity = eps
In [49]:
rce.integrate_years(2.)
In [50]:
# Check for energy balance
rce.diagnostics['ASR'] - rce.diagnostics['OLR']
Out[50]:
In [51]:
plot_sounding([col, rce])
Out[51]:
In [52]:
# ANother 1% increase in absorptivity
rce2 = climlab.process_like(rce)
rce2.subprocess['LW'].absorptivity *= 1.02
In [53]:
rce2.compute_diagnostics()
RF = -(rce2.diagnostics['OLR'] - rce.diagnostics['OLR'])
print 'The radiative forcing is %f W/m2.' %RF
In [54]:
# Timestep forward, and the check for energy balance
rce2.integrate_years(2.)
rce2.diagnostics['ASR'] - rce2.diagnostics['OLR']
Out[54]:
In [55]:
plot_sounding([col, rce, rce2])
Out[55]:
In [56]:
ECS = rce2.Ts - rce.Ts
print 'Equilibrium climate sensitivity is %f K.' %ECS
In [57]:
# Calculate the net climate feedback
# This is the change in TOA flux per degree warming that was necessary to get back to equilibrium.
feedback = -RF/ECS
print 'The net feedback is %f W/m2/K' %feedback
In [58]:
# could calculate a Planck feedback explicitly...
# What would the TOA flux change be if the warming were perfectly uniform?
rce3 = climlab.process_like(rce)
rce3.subprocess['LW'].absorptivity *= 1.02
rce3.Ts += ECS
rce3.Tatm += ECS